function F = cal_firm_vars(param, const, eqm)

    % Useful function    
    firmvar = @(var) var(eqm.q_index)';
    
    
%% Quality, Productivity

    % Quality level
    firm.qlevel = firmvar(const.Qgrid);
    
    % Productivity
    firm.lnz = full(sum(const.lnzQ.*eqm.q_indicator, 2));

    
%% Export probabilities
    
    firm.prob_exp = full(sum(eqm.ExportProbQ.*eqm.q_indicator, 2));    
    firm.prob_nonexp = full(sum((1-eqm.ExportProbQ).*eqm.q_indicator, 2));
    
    
%% Export intensity    
    
   
    firm.exp_intst_exp = firmvar(1-eqm.rv_rev);

        
    
%% Sales, Indegree, Outdegree

    % Log sales    
    firm.lnsales_exp = (param.sigma-1)*const.gamma*firm.lnz + log(firmvar(eqm.Pi1));    
    firm.lnsales_nonexp = (param.sigma-1)*const.gamma*firm.lnz + log(firmvar(eqm.Pi0)); 
        
    % Number of suppliers
    firm.indegree_exp = exp(firm.lnsales_exp).^(1/param.beta_m).*firmvar(const.C_m.*eqm.theta_m);
    firm.indegree_nonexp = exp(firm.lnsales_nonexp).^(1/param.beta_m).*firmvar(const.C_m.*eqm.theta_m);
       
    % Number of customers
    tempsum = sum(const.phi_v.*kron(ones(1,param.Q),eqm.theta_v'), 1);
    firm.outdegree_exp = exp(firm.lnsales_exp).^(1/param.beta_v).*firmvar(eqm.rv_ads.*eqm.C_v1.*tempsum);
    firm.outdegree_nonexp = exp(firm.lnsales_nonexp).^(1/param.beta_v).*firmvar(const.C_v0.*tempsum);
    
    
%% Wage, Employment

    % Average wage per worker
    firm.wage = firmvar(param.w.*const.wage_grid);    

    % Employment (number of workers)
    frac = (1-param.alpha_m-param.alpha_s)*(param.sigma-1)/param.sigma;
    firm.employment_exp = frac.*exp(firm.lnsales_exp)./firm.wage;
    firm.employment_nonexp = frac.*exp(firm.lnsales_nonexp)./firm.wage;      
    
    
%% Average log wage of suppliers

    % Unweighted average supplier wage
    avg_unwgt_lnwageS_q = 1./eqm.V.*sum(const.phi_v.*kron(ones(param.Q,1), ...
        (log(param.w)+const.lnwage_grid).*(const.C_v0.*eqm.Pi0.^(1/param.beta_v).*eqm.integ.E_zv0 ...
        +eqm.rv_ads.*eqm.C_v1.*eqm.Pi1.^(1/param.beta_v).*eqm.integ.E_zv1)),2)';
    firm.avg_unwgt_lnwageS = firmvar(avg_unwgt_lnwageS_q);
    
    % Expenditure weighted average supplier wage
    avg_wgt_lnwageS_q = eqm.theta_m.*eqm.c.^(param.sigma-1)./eqm.V.*sum(const.phi_y...
        .*const.phi_v.*kron(ones(param.Q,1), (log(param.w)+const.lnwage_grid).*eqm.Psig),2)';    
    firm.avg_wgt_lnwageS = firmvar(avg_wgt_lnwageS_q);
      
    
%% Fraction of higher quality match

    % Fraction of higher (or same) quality supplier
    flag_S_higher_q = zeros(param.Q, param.Q); %row is buyer, column is supplier
    for i = 1:param.Q
        flag_S_higher_q(i,:) = [zeros(1,i-1), ones(1,param.Q-i+1)];
    end
    frac_S_higher_q = sum(flag_S_higher_q.*const.phi_v.*repmat(eqm.Vbar,param.Q,1),2)'./eqm.V;
    firm.frac_S_higher_q = firmvar(frac_S_higher_q);
    
    % Fraction of higher quality customer
    flag_B_higher_q = zeros(param.Q, param.Q); %row is buyer, column is seller
    for i = 1:param.Q
        flag_B_higher_q(:,i) = [zeros(i-1,1); ones(param.Q-i+1,1)];
    end    
    tempmatrix = const.phi_v.*repmat(eqm.theta_v',1,param.Q);
    frac_B_higher_q = sum(flag_B_higher_q.*tempmatrix,1)./sum(tempmatrix,1);        
    firm.frac_B_higher_q = firmvar(frac_B_higher_q);
            

%% Return

    F = firm; 

    
end